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We derive simple analytical expressions for the particle 
density p(r) and the kinetic energy density r(r) for a system of 
noninteracting fermions in a d-dimensional isotropic harmonic 
oscillator potential. We test the Thomas-Fermi (TF, or local- 
density) approximation for the functional relation r[p] using 
the exact p(r) and show that it locally reproduces the exact 
kinetic energy density r(r), including the shell oscillations, 
surprisingly well everywhere except near the classical turning 
point. For the special case of two dimensions (2D), we obtain 
the unexpected analytical result that the integral of ttf [p(r)] 
yields the exact total kinetic energy. 

PACS numbers: 03.65.Sq, 05.30.Fk, 31.15.Ew, 71.10.Ca 

The evaporative cooling of dilute (i.e., almost nonin- 
teracting) fermionic gases has recently been achieved by 
Jin and DeMarco at JILA in Colorado This spec- 
tacular experimental milestone has stimulated an enor- 
mous effort to explore and understand the properties of 
these new quantum systems, which can be viewed as 
the quantum analogue of the Bose-Einstein condensation 
(BEC) recently observed in ultracold trapped Bose gases. 
While it is true that the theory of homogeneous dilute 
Fermi systems is fairly well developed, addressing new 
experiments probing strongly inhomogeneous systems in 
regimes far from equilibrium, represents a much greater 
challenge to theory. To this end, Vignolo et al. Q have 
recently used a Green's function method to compute the 
particle and kinetic energy densities for a system of non- 
interacting fermions in a one-dimensional (ID) harmonic 
oscillator potential. In principle such a (quasi)-lD sys- 
tem can be achieved experimentally using state-of-the- 
art magnetic confinement techniques jjj. Owing to the 
enhanced shell structure found in ID, Vignolo et al. sug- 
gest that these quantum oscillations may be accessible to 
observation in magnetically trapped gases of fermionic al- 
kali atoms. We give here much simpler analytical results 
for the more general case of a d-dimcnsional harmonic 
potential and use them to test their Thomas-Fermi (TF) 
functional relation. 

The work presented in this paper is also applicable to 
a 2D electron gas confined to so-called quantum dots || . 
The external confinement potential of these dots is in 
many cases essentially harmonic. Bhaduri et al. Q have 
shown that the inclusion of a short-range two-body inter- 
action may be included via fractional statistics, provided 



that one uses the TF relation ttf[/o] relevant for 2D. 

The method. — We start from a system of noninter- 
acting fermions described by the time-independent 
Schrodinger equation 



H</>i{r)= f + V(r) 0i(r) = ei &(r) 



(1) 



where V(r) is a local potential to be specified later. The 
single-particle density matrix can be obtained by an in- 
verse Laplace transform of the Bloch density matrix: 



p(r,r') = 2]T tf(r')&(r)=£ii 

6i<E F 



, (2) 



where the latter is defined by 

C(r,r';/3) = £#(r')&(r) exp{-/3ej. (3) 

all i 

Ep is the Fermi energy, the factor 2 accounts for spin. 
We now use center-of-mass and relative coordinates: 



q=-(r + r'), s = r - r', 



(4) 



so that the local density is p(q) = p(q, s)| s=0 - For the 
kinetic energy density, we investigate two expressions || : 



r(q) 



2m 



2^ #(q)VVi(q), 



ri (q) = — 2£ |V^(q)| 2 
2m ^— ' 



(5) 
(6) 



In the presence of time-reversal symmetry they are sim- 
ply related by 



A convenient quantity is their mean, 



1 



C(q) = 2[r(q)+n(q)], 
which is obtained from the density matrix by 



£(q) 



2m 



[V^(q )S )] s=0 , 



(7) 



(8) 



(9) 



1 



where V s is the gradient with respect to the variable 
s. Note that all three quantities r(q), n(q), and £(q) 
integrate to the same exact kinetic energy. 

We now specialize to an isotropic harmonic oscillator 
potential in d dimensions: 



V(r) 



(10) 



where r = \J x\ + . . 
exact Bloch density matrix for this system is given by 



x 2 is the radial variable. The 



C(r,r';p)=C(q,s;p) 



exp 



raw 



q tanh 



'f3huj 



d/2 



1 



smh d/2 {[3hu) 



-ctgh 



(11) 



To get the particle and kinetic energy densities, we need 
to perform inverse Laplace transforms of the above func- 
tion and its derivatives at s = 0. For the first exponential 
factor in (|ll]) , we employ the following relation which can 
be derived from Ref. 

exp{-xtanh(/3/2)} 



E 

n=0 



(-l) n L n (2x) e~ x {e- nP + e -(»+ 1 )/ 3 | . (12) 



This relation holds if \z\ — \e~^\ < 1, which is fulfilled 
since the contour of the inverse Laplace transform inte- 
gral in the complex (3 plane goes along j3 = it + c for 
t € (— oo, oo) with c > 0. Further analytical progress 
depends on the dimensionality d. 

The case d — 2. — Here we can directly use the fol- 
lowing exact Laplace inverse M 



-n/3 



2 V 6[A- (2k 



k=0 



1) 



(13) 



_/3sinh(/3)_ 

When filling M + 1 oscillator shells, the Fermi energy is 

E F = Hlo (M + 1 + 5), (14) 

with S being an infinitcsimally small positive number. 
Combining Eqs. (|l^), ( |l3| ) and carefully evaluating the 
sums over the step functions, we get for the density 



M 

p(q) = 2 (^) ^(M-/x+l)(-l)^(2x) e 

p,=0 



(15) 



where x = (niLo/Ti)q 2 . The kinetic energy density (||) 
is given, after some suitable manipulations of hyperbolic 
functions, by the following Laplace inverse: 



(muj 
huj I — 

V TTh, 



(3 2sinh J ((3hw/2) 



exp {— xta.nh(f3Tiuj/2)} 



(16) 



Removing one inverse sinh factor by the identity 



2swh(f3hu/2 



1 OO 

i = v P -( m+1 



(17) 



and proceeding as above, we get the final expression for 
the kinetic energy density: 

M 

- huj (^) Y^tM-^ + lfi-lYL^x) e~\ (18) 

The integrals d 2 q of the densities ([l5|), ( |lS| ) are readily 
evaluated using 



L n (2x)e~ x dx = (-1)", 



(19) 



and yield the correct results for the number N of particles 
in M + 1 filled shells 

„ M 

p(q)d 2 q = 2^2^+1) = M 2 + 3M + 2 = N(M) , (20) 

and for their exact kinetic energy E^in (M) 
„ M 

/ e(«) a = 7k^( m +i) 2 



() 



(2M 3 + 9M 2 + 13M + 6) = E kin (M) . (21) 



Next we investigate the Thomas-Fermi (TF) relation 
between r (or t\ or £) and p, which in 2D is 



t tf[p] = —irp 2 , 
2m 



(22) 



see also Eq. fl33| ) below. Inserting Eq. ( |15| ) into the right- 
hand side above and integrating, using the orthonormal- 
ity of the Laguerre polynomials, we find 

„ M 

/ T TF [p{q)] d 2 q = + ^ = E Mn(M) . (23) 

This means that the simple TF functional - without gra- 
dient corrections Q - using the exact density p(q) yields 
the exact quantum-mechanical kinetic energy, which is 
highly nontrivial and unexpected. The local behavior of 
T TF[p{q)\ will be examined numerically below. 

The cases d — 1 and d > 3. — For d = 1, we have a 
square root of the sinh factor in the denominator of (|ll]) . 
To handle it, we use the expansion 



sinh 1/2 (s) 



V2 



, ) 1/J = -5=e'/Vl-e- J « 



v ^ rn—1 



(2m -3)!! _ 5 
(2m)!! 6 



. (24) 



2 



which converges for Re s > 0, and include it as a factor 
on top of the d = 2 case. For d = 3, we have to include 
its inverse, which has the expansion 



sinh- 1 / 2 ( s ) = V2e- S / 2 (1- 



For d = 4, we need 

oo 

sinh-^s) = 2(e s - e^) -1 = 2 J] 



-(2m+l)s 



(25) 



(26) 



m= 



and so on. Using Ep = hto (M + d/2 + 5) and proceeding 
as above for the Laplace inversions, we obtain the general 
expressions 



ti=0 



(27) 
(28) 



The coefficients Fi rf) , are given by 



[.z/2] 



f + iy + 1 - 2m ) 



[u/2] 



G« = [v + l) 2 + £ (i/ + f - 2m) 2 . 9 W , (29) 

m— 1 

using (for m > 1) 

g g) = -(2m - 3)!!/(2m)!! with .9^ = 1/2, 
£ 2) =0, </£> = (2m-l)H/(2ro)!!, 



-,(4) 



1 



-,(5) - 



(2m + l)!!/(2m)H, 



(30) 



and so on. For even dimensions, Fjfi and Gu can be 
computed analytically, with [v/2] = integer (i//2). 

Like Eqs. @, @, to which Eqs. @, @ reduce for 
2D, the latter are much simpler for numerical computa- 
tions than their definitions in terms of the eigenfunctions, 
which necessitate multiple summations for d > 2. 

In the TF or local-density approximation (LDA), one 
gets from the potential V(q) = (mui 2 /2)q 2 — (hoj/2)x 
the following local densities 



4 f fmuo\ d / 2 



dr(|) 



(31) 



= (5+2) r(l) (2^) ( A --/ 2 ) d/2+1 ' ( 32 ) 

where A = Ep/(hu>). The TF functional relation between 
r and p is then given by 



t T f [p] 



h z And 
2m (d + 2) 



dU 
4 \2 



2/d 



.1+2/ d 



(33) 



Our objective now is to study numerically the above re- 
lation using the exact densities and to see how well it 
holds locally as well as globally (i.e., upon integration). 
We have already seen analytically that its integral yields 
the exact kinetic energy for d = 2. 

Numerical results. — The following figures show nu- 
merical results in units such that % = lu = m = 1. We 
make the following observations: 

1. As is well known (cf. also Ref. Q), the densities p(q) 
and r(q), T\(q) oscillate around the smooth TF densities 
( |3l| ) and j32|) , respectively, except near the turning point 
where the latter go to zero (see Figs. ||, ||). 

2. The shell oscillations in the quantities r(q) and T\ (q) 
are exactly opposite in phase, so that their mean £(q) is a 
smooth function of q that (except near the turning point) 
closely follows the TF density TTF^q) given in ( |32|) (see 
Fig. [I]). This has already been observed long ago 

3. The functional ttf [/>(</)] in (^3|), using the exact 
densities p(q), reproduces locally the exact kinetic energy 
density r(q) surprisingly well, including the shell oscilla- 
tions, except near the classical turning point and in the 
far tail region (see Figs. ^, ||). 

4. The integral of t>t.f (/>(<?)] over the d-dimensional 
space yields the exact kinetic energy only for d = 2, see 
Eq. (pi) . In the other cases it yields kinetic energies with 
errors less than one percent for N > 14, 100, and 900 in 
d = 1, 3, and 4 dimensions, respectively. 
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FIG. 1. Kinetic energy densities for N = 240 particles 
filling 8 shells of a 3D isotropic harmonic oscillator. Upper 
panel: solid line: r(q), dashed line: ri(q), dotted line: £((/)• 
Lower panel: solid line: £(<?), dashed line: TTF(q). 
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FIG. 2. Kinetic energy densities for N = 50 particles fil- 
ling 25 shells of a ID harmonic oscillator. Solid lines: exact 

using the exact 
The inserts give 



r(q). Dashed lines: TF relation ttf \p(q)1 
p(q). Dotted lines: TF density ttf(q) 
close-ups near the center and the tail region. 



reproduces the strong local shell oscillations in r(q) so 
accurately - in the figures, the error cannot been rec- 
ognized except in the tail regions - is therefore unex- 
pected and does not seem to have been noticed before 
plj. We must, however, add the caveat that the func- 
tional ttf[p(q)] cannot be variationally exact. As is well 
known, indeed, the Euler-Lagrange variational equation 
derived from it leads precisely to the TF density ptf{q) 
in Eq. and not to the exact quantum-mechanical 

density p(q). 

In the 2D case, where the integral even reproduces the 
exact kinetic energy, our result supports the basic as- 
sumptions made in Ref. (J] concerning the inclusion of a 
short-range two-body force through fractional statistics, 
which relies upon the TF relation (0). 

Finally, we wish to emphasize that the recent work of 
Vignolo et al. || is a special case of our more general 
results, and point out that the prominent shell structure 
displayed in 2D could also become observable in experi- 
ments on alkali vapours. 
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FIG. 3. The same as Fig. |, but for TV = 420 particles 
filling 20 shells of a 2D isotropic harmonic oscillator. 

Summary and conclusions. - Our finding that the 
TF functional relation Ttf[p(<i)] works so well locally 
is surprising, since it theoretically is exact only in the 
LDA, i.e., for spatially homogeneous systems. That it 
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